module polynom
use technical
use df_dt
implicit none
public 
contains
!*********************************************
subroutine poly(ya,xa,dim,x,y,dy)
! polynomial extrapolation
real, dimension(dim) :: ya, xa
real :: x,y,dy
integer :: dim
real :: den, dif, dift, ho, hp, w
real, dimension(dim) :: c, d
integer :: i,m,ns=1
!-----------------------------------------
dif=abs(x-xa(1))
do i=1,dim
   dift=abs(x-xa(1))
   if (dift < dif) then
      ns=i
      dif=dift
   endif
enddo
c=ya
d=ya
y=ya(ns)
do m=1,dim
   do i=1,dim-m
      ho=xa(i)-x
      hp=xa(i+m)-x
      w=c(i+1)-d(i)
      den=ho-hp
      if (den==0) call fatal_error('subroutine poly','Two points identical')
      den=w/den
      d(i)=hp*den
      c(i)=ho*den
   enddo
   if (2*ns < (dim-m)) then
      dy=c(ns+1)
   else 
      dy=d(ns)
   endif
   y=y+dy
enddo
end subroutine poly
!********************************
subroutine stoer(y,dydx,nv,xx,htry, eps, yscal, hdid,hnext)
! Bulirsch-Stoer extrapolation
implicit none
integer, parameter :: kmaxx=8
real, parameter :: safe1=0.25, safe2=0.7,scalmx=0.1
real, parameter :: redmax=1.0e-5, redmin=0.7, tiny=1.0e-30
integer, parameter :: imaxx=kmaxx+1 
!
integer :: nv
integer :: i,iq,k,kk,km
integer :: first=1,kmax,kopt
real, dimension(nv) :: y, dydx, yscal
real :: xx, htry,eps,hdid
real ::  epsold = -1.0,xnew=-1.0e29
real:: eps1,errmax,fact,h,red,scale,work,wrkmin,xest,hnext
real, dimension(kmaxx) :: err,x
real, dimension(imaxx+1) :: a
real, dimension(KMAXX+1,KMAXX+1) :: alf
real, dimension(nv,kmaxx) :: d
real, dimension(nv) :: yerr,ysav,yseq
integer, dimension(IMAXX+1) :: nseq =(/0,2,4,6,8,10,12,14,16,18/)
integer ::  reduct,exitflag=0;

!
!----------------------------------
if (eps == epsold) then
   hnext = xnew
   eps1=SAFE1*eps
   a(1)=nseq(1)+1
   do k=1,KMAXX  
      a(k+1)=a(k)+nseq(k+1)
   enddo
   do iq=2,KMAXX
      do k=1,iq-1
         alf(k,iq)=eps1**((a(k+1)-a(iq+1))/((a(iq+1)-a(1)+1.0)*(2*k+1)))
      enddo
   enddo
   epsold=eps
   do kopt=2,KMAXX-1 
      if (a(kopt+1) > a(kopt)*alf(kopt-1,kopt)) exit
   enddo
   kmax=kopt
endif
h=htry
ysav=y
if (xx /= xnew .or. h /= hnext) then 
   first=1
   kopt=kmax
endif
reduct=0
do
   do k=1,kmax
      xnew=xx+h
      if (xnew == xx) call fatal_error('subroutine stoer','step size underflow')
      call mmid(ysav,dydx,nv,xx,h,nseq(k),yseq,x,d,kmaxx)
      xest=(h/nseq(k))**2
      call pzetra(k,xest,yseq,y,yerr,nv,x,d,kmaxx)
      if (k /= 1) then
         errmax=tiny
         do i=1,nv 
            errmax=max(errmax,abs(yerr(i)/yscal(i)))
         enddo
         errmax=errmax/eps
         km=k-1
         err(km)=(errmax/safe1)**(1.0/(2*km+1))
      endif
      if (k/=1 .and. ((k>=kopt-1) .or. (k>=first))) then
         if (errmax <1.0) then
            exitflag=1
            exit
         endif
         if ((k==kmax) .or. (k==kopt+1)) then
            red=safe2/err(km)
            exit
         else if ((k==kopt) .and. (alf(kopt-1,kopt)< err(km))) then
            red=1.0/err(km)
            exit
         else if ((kopt==kmax) .and. (alf(km,kmax-1) < err(km))) then
            red=alf(km,kopt)/err(km)
            exit
         else if (alf(km,kopt) < err(km)) then
            red=alf(km,kopt-1)/err(km)
            exit
         endif
      endif
   enddo
   if (exitflag/=0) exit
   red=min(red,redmin)
   red=max(red,redmax)
   h=h*red
   reduct=1
enddo
xx=xnew
hdid=h
wrkmin=1.0e35
do kk=1,km
   fact=max(err(kk),scalmx)
   work=fact*a(kk+1)
   if (work<wrkmin) then
      scale=fact
      wrkmin=work
      kopt=kk+1
   endif
enddo
hnext=h/scale
if ((kopt >=k) .and. (kopt /= kmax) .and. reduct==0) then
   fact=max(scale/alf(kopt-1,kopt),scalmx)
   if (a(kopt+1)*fact <=wrkmin) then
      hnext=h/fact
      kopt=kopt+1
   endif
endif
end subroutine stoer
!*******************************************
subroutine mmid(y,dydx,nvar,xs,htot,nstep,yout,x,d,kmaxx)
! midpoint integration
implicit none
Real, dimension(kmaxx) :: x
real, dimension(nvar,kmaxx) :: d
real, dimension(nvar) :: y, dydx,yout
real :: xs,htot,h,h2,swap
integer :: nvar,nstep,kmaxx
integer :: n,i
real, dimension(nvar) :: ym,yn
!-------------------------------------
h=htot/nstep
ym=y
yn=y+h*dydx
x=xs+h
call eval_rhs(x,yn,yout)
h2=2.0*h
do n=2,nstep
   do i=1,nvar
      swap=ym(i)+h2*yout(i)
      ym(i)=yn(i)
      yn(i)=swap
   enddo
   x=x+h
   call eval_rhs(x,yn,yout);
enddo
yout=0.5*(ym+yn+h*yout)
end subroutine mmid
!************************************************
subroutine pzetra(iest,xest,yest,yz,dy,nv,x,d,kmaxx)
implicit none
real, dimension(nv) :: yest,yz,dy,c
real, dimension(kmaxx) :: x
real, dimension(nv,kmaxx) :: d
real :: xest,q,f2,f1,delta
integer :: iest,nv, k,hj,kmaxx,j,k1
!----------------------------
x(iest)=xest
dy=yest
yz=yest
if (iest==1) then
   do j=1,nv
      d(j,1)=yest(j)
   enddo
else 
   c=yest
   do k1=1,iest-1
      delta=1.0/(x(iest-k1)-xest)
      f1=xest*delta
      f2=x(iest-k1)*delta
      do j=1,nv
         q=d(j,1)
         d(j,k1)=dy(j)
         delta=c(j)-q
         dy(j)=f1*delta
         c(j)=f2*delta
         yz(j)=yz(j)+dy(j)
      enddo
   enddo
endif
end subroutine pzetra
!
end module

